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The experimental two-microphone transfer function method (TMTF) is adapted to the numerical 
framework to compute the radiation and input impedances of three-dimensional vocal tracts of 
elliptical cross section. In its simplest version, the TMTF method only requires measuring the 
acoustic pressure at two points in an impedance duct and the postprocessing of the corresponding 
transfer function. However, some considerations are to be taken into account when using the TMTF 
method in the numerical context, which constitute the main objective of this paper. In particular, 
the importance of including absorption at the impedance duct walls to avoid lengthy numerical 
simulations is discussed and analytical complex axial wave numbers for elliptical ducts are derived 
for this purpose. It is also shown how the plane wave restriction of the TMTF method can be 
circumvented to some extent by appropriate location of the virtual microphones, thus extending 
the method frequency range of validity. Virtual microphone spacing is also discussed on the basis 
of the so called singularity factor. Numerical examples include the computation of the radiation 
impedance of vowels /a/, /i/ and /u/ and the input impedance of vowel /a/, for simplified vocal 
tracts of circular and elliptical cross sections. 
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I. INTRODUCTION 

Radiation and input impedances of vocal tracts are mag- 
nitudes of special interest for voice production. They 
play a significant role in determining wave radiation from 
the mouth or in modeling the acoustic coupling between 
vocal tract and vocal folds acoustics. Impedances can 
be computed from numerical simulations of vocal tract 
acoustics. Due to vocal tract intricate geometry, the 
most extended numerical approach to carry out these 
simulations is the finite element method (FEM). Several 
works can be found in literature dealing with FEM com- 
putations both in the frequency domain (e.g., Matsuzaki 
et al, 2000; Motoki, 2002; Hannukainen et a/., 2007) and 
time domain(e.g., Svancara and Horacek, 2006; Vampola 
et al, 2008, 2011). Ocasionally other approaches such 
as finite differences have also been used (e.g., Takemoto 
et al, 2010). 

With regards to voice production, the natural choice 
turns to be that of working in the time domain. If the 
wave equation is solved in its mixed form (e.g., Takemoto 
et al, 2010; Codina, 2008), impedances can be directly 
computed from the Fourier transforms of the acoustic 
pressure and acoustic velocity time evolutions. However, 
this is not possible if the wave equation is solved in ir- 
reducible form for the acoustic pressure (e.g. Vampola 
et a/., 2011) or for the velocity potential (e.g., Matsuzaki 
et a/., 2000). In such cases the acoustic velocity has to 
be computed from the acoustic pressure or the velocity 
potential gradients, the convergence error being of higher 
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order than for the primary variables (e.g., Hughes, 2000). 

In this paper an alternative is proposed to compute vo- 
cal tract impedances from time domain simulations with- 
out having to compute the acoustic velocity field. The 
idea is to adapt the experimental two-microphone trans- 
fer function method (TMTF) to the numerical frame- 
work. Originally developed by Chung and Blaser (1980), 
the simplest version of the experimental TMTF method 
only requires measuring the time evolution of the acoustic 
pressure at two points in an impedance duct, and com- 
puting the corresponding transfer function. From this 
transfer function, the radiation and/or input impedances 
at a given surface can be derived. 

Although at first sight applying the TMTF method 
to compute vocal tract impedances may look straightfor- 
ward, some issues are worth exploring. These constitute 
the basis of this work. On the one hand, we will show 
the importance of dealing with lossy impedance ducts to 
reduce the overall time duration of the simulations. The 
inclusion of wall losses allows to strongly attenuate the 
duct first eigenmode, which otherwise determines the to- 
tal duration of the computation. However, this implies 
using appropriate complex wavenumbers in the TMTF 
expressions, which are well-known for three-dimensional 
circular cylindrical ducts, and which will be derived in 
this work for elliptical cross sectional impedance ducts, 
given their importance in voice production (see e.g., Mo- 
toki, 2002; Matsuzaki et al, 2000, where elliptical vocal 
tracts are used). On the other hand, the frequency range 
of validity of the TMTF method will be analyzed. First, 
it will be shown how the plane wave propagation restric- 
tion that limits the TMTF maximum frequency, can be 
generously surpassed by appropriate location of the vir- 
tual microphones (mesh nodes where the acoustic pres- 
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sure time evolution is collected). Second, the appropriate 
virtual microphone spacing will be determined by means 
of the so called singularity factor (SF) introduced by Jang 
and Ih (1998), according to the maximum frequency of 
analysis. Throughout the work, time domain FEM sim- 
ulations for the irreducible wave equation will be per- 
formed with an in- house software, to compute vocal tract 
impedances using the adapted TMTF method. However, 
any other time domain numerical approach could benefit 
from the hereafter exposed results. A preliminary version 
of some of them was presented by the authors in Arnela 
and Guasch (2012). 

The paper is organized as follows. Section II presents 
the methodology followed to compute the acoustic 
impedance of vocal tracts from numerical simulations. 
In section III, the various considerations to be taken into 
account when adapting the TMTF method to the numer- 
ical framework become analyzed. Numerical examples of 
computed impedances for vocal tracts of circular and el- 
liptical cross section are provided in section IV. Finally, 
conclusions close the paper in section V. 



II. METHODOLOGY 

A. The two-microphone transfer function method 

A brief survey of the two-microphone transfer function 
method (TMTF) will be next presented. This method 
was originally developed by Chung and Blaser (1980) and 
later on standardized in the ISO 10534-2 (1998) for mea- 
suring the normal reflection coefficient of material sam- 
ples. The specific acoustic impedance Z can be obtained 
from the latter. The TMTF proceeds as follows. First, 
plane waves are generated at the entrance of a duct of 
length L, referred to as the impedance duct. The acous- 
tic pressure signals P\(f) and P2(f) become measured 
at two points x\ and x 2 close to the impedance duct 
exit, which is designated as the reference surface, and 
the transfer function H\2(f) = P2(f)/Pi(f) is computed. 
The reflection coefficient R\ at position x\ is given by 



Ri = 



H 



12 



Hr — H] 



12 



(i) 



with Hi and Hr respectively standing for the incident 
and reflected wave transfer functions. Assuming plane 
wave propagation, the transfer functions Hj and Hr be- 
come Hi = e~ jkzS , H R = e jkzS , with s = \x± — x 2 \ being 
the distance between microphones and k z the wavenum- 
ber in the axial direction. In order to translate the reflec- 
tion coefficient to the reference surface, defined at x = 0, 
a factor e^ 2kzXl is introduced in (1), leading to the fol- 
lowing expression for the normal reflection coefficient 



R = R ie j2kzXl = 
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(2) 



The normalized specific acoustic impedance Z' can be 
finally obtained by means of 



Z' = Z/Z = 



R 



(3) 
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FIG. 1. (Color online) A sketch for problem (4) in text. Com- 
putational domain for the input impedance (left) and for the 
radiation impedance (right). 



with Z standing for the specific acoustic impedance and 
Zq = PqCq for the characteristic impedance. Z is usually 
split in its real R (resistive) and imaginary X (reactive) 
components, Z = Zq(R + jX). 



B. Problem statement 

Suppose that we want to compute the radiation 
impedance Z r or the input impedance Z in of a vocal 
tract from a time domain numerical simulation using 
the two-microphone transfer function (TMTF) method. 
We first need to couple an impedance duct to the ref- 
erence surface, the duct having the same section and 
shape than this surface. For instance, to compute the 
input impedance of a vocal tract, the impedance duct 
will be coupled to the vocal tract at the glottal section 
(see Fig. 1), whilst for the radiation impedance, the duct 
will be coupled to the mouth exit replacing the existing 
vocal tract geometry (see Fig. 1). 

The next step consists in carrying out a time do- 
main numerical simulation for the acoustic pressure evo- 
lution. Let us denote by Q the finite computational do- 
main with boundary dQ. We can identify on dft four 
non-intersecting regions (see Fig. 1): Tq and Tz for the 
impedance duct boundaries, Tw for the vocal tract walls, 
Th for the human head boundary, and for the exter- 
nal boundary, where a fictitious non-reflecting condition 
has to be imposed. We aim at finding the acoustic pres- 
sure field p(x,t) in Q that satisfies 







(d 2 tt -c 2 V 2 )p{x,t) 
with boundary and initial conditions 



in ft, t > 0, 
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in t 


= 0. 


(4g) 



In Eq. (4), Co denotes the speed of sound, po the air 
density, S the impedance duct cross sectional area at Tq, 
dt = d/dt designates the partial time derivative and n 



M. Arnela and O. Guasch.: Elliptical vocal tract impedance computation 2 



FIG. 2. (Color online) Geometries used for computing the in- 
put impedance (right top) and the radiation impedance (right 
bottom) of the elliptical vowel /a/ (left). 



FIG. 3. (Color online) Cut with surface mesh details of the 
impedance duct for the elliptical /a/ radiation impedance 
computation. Dots indicate locations of virtual microphones 
for capturing Pi(f) and P2(f) in text. 



the normal vector pointing outwards a surface. With re- 
gards to boundary conditions (BC), Q(t) in (4b) stands 
for a volume velocity generated by an imaginary loud- 
speaker. The BCs in (4c) and (4d) account for constant 
frequency losses at the inner walls, being /i the boundary 
admittance coefficient (subindexes w and z simply indi- 
cate that different absorption values can be introduced 
at each boundary). \i is related to the wall impedance 
Z w by means of /i = pqCq/Z w . The BC (4e) expresses 
that the human head is taken as a rigid surface (/i = 0). 
Finally, BC (4f) is the well-known Sommerfeld radiation 
condition, which guarantees that emanating waves from 
the mouth propagate outwards to infinity. However, this 
condition is only optimal for sound waves impinging on 
in the normal direction. To overcome this problem, 
we have made use of a Perfectly Matched Layer (PML) 
(Berenger, 1994) for the wave equation in its irreducible 
form. In particular, we have adapted the PML origi- 
nally developed for the finite difference framework (Grote 
and Sim, 2010) and formulate it for our finite element 
code. Details on the implemented numerical scheme can 
be found in appendix A. 

As the simulation evolves, the acoustic pressure sig- 
nals pi(t) and P2(t) are collected at two nodes of the 
finite element mesh, as if they were virtual microphones. 
From their Fourier Transform we obtain Pi (/) and P2(/)> 
and making use of the TMTF method described in sec- 
tion II. A, we can compute the acoustic impedances. 

C. Vocal tract models 

The vocal tracts of the three extreme cardinal vowels 
/a/, /i/ and /u/ with circular and elliptical cross section 
have been considered. In order to shorten notation, we 
will refer to them for instance as circular /a/ or ellipti- 
cal /a/. For the circular vocal tracts, we have used the 
simplified vocal tract geometries generated from sections 
provided by Story (2008). With regards to the elliptical 
vocal tracts (see Fig. 2), we have reshaped the circu- 
lar sections according to the eccentricity of the elliptical 
mouth apertures described in Fromkin (1964). A spheri- 
cal surface of radius 0.09 m has been used to emulate the 



human head in all cases. 

Following the procedure outlined at the beginning 
of section ILB (see Fig. 2), we have replaced the vo- 
cal tract geometry by an impedance duct with equal 
mouth aperture to compute the vocal tract radiation 
impedance, whereas we have coupled an impedance duct 
at the glottal section of the vocal tract to compute its 
input impedance. The impedance duct has a length of 
L = 0.1 m to fulfill the requirements of the standard 
ISO 10534-2 (1998) (the length should be at least three 
times the duct radius or the major semi-axis). The vir- 
tual microphones have been located at the centerline of 
the impedance duct and separated a distance s = 0.01 m 
apart (see Fig. 3). Following the recommendations of the 
ISO 10534-2 (1998), the first virtual microphone has been 
placed at a distance from the reference surface slightly 
larger than two times the impedance duct radius, or the 
major semi-axis. Concerning the reference surface for vo- 
cal tract radiation impedances computations, it should 
be noted that it is well-defined for the circular case given 
that the intersection of a cylindrical vocal tract with a 
spherical human head results in a flat disk. However, this 
is not the case if an elliptical vocal tract is used. For such 
cases we have chosen the elliptical section where the ma- 
jor semi- axis intersects the sphere as the reference surface 
(see Fig. 3). 

D. Simulation details 

The computational domain consists of an outer volume 
of dimensions 0.25 x 0.2 x 0.2 m, where the spherical 
head has been placed so that sound waves can emanate 
from the mouth. This volume has been surrounded with 
a PML of width 0.1 m to absorb any incident wave. The 
PML has been configured to get a relative reflection co- 
efficient of ^ = 10 -4 (see appendix A). The computa- 
tional domain has been meshed using unstructured tetra- 
hedral elements with a size comprising h « 0.1 cm inside 
the impedance duct, h ~ 0.5 cm in the outer volume and 
h w 0.75 cm in the PML region (see Fig. 3 to appreciate 
some mesh details). 
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Problem (4) with a PML included has been solved us- 
ing the finite element approach described in appendix A. 
A time interval of total duration T = 30 ms has been sim- 
ulated using a sampling rate of f 3 = 1/At = 2000 KHz. 
The values Co = 345 m/s and po = 1.1933 kg/m 3 have 
respectively been chosen for the speed of sound and for 
the air density. Concerning boundary conditions, a wide- 
band impulse has been used for the volume velocity Q(t) 
in (4b), and input at the impedance duct entrance. We 
have used a gaussian pulse of the type (Takemoto et al., 
2010) 



gp(n) 



,[(At 



)0.29T q , 



[m 3 /s], 



(5) 



with T gp = 0.646//o and fo = 10 KHz. To avoid nu- 
merical errors beyond the maximum frequency of inter- 
est (fmax — 10 KHz), this pulse has been filtered using 
a low-pass filter with cutoff frequency 10 KHz. For the 
boundary admittance coefficient at the vocal tract walls 
we have used \i w = 0.005, which corresponds to the wall 
impedance of the vocal tract tissue Z w = 83666 Kg/m 2 s 
(see Svancara and Horacek, 2006). For the impedance 
duct, we have used the artificial value ji z =0.01. 



III. THE TWO-MICROPHONE TRANSFER FUNCTION 
METHOD FOR NUMERICAL SIMULATIONS 

A. Damping the first impedance duct eigenmode 

Although theoretically numerical simulations to compute 
vocal tract impedances could be carried out using a loss- 
less impedance duct, including boundary losses is manda- 
tory for the simulations to have a reasonable duration. 
From an experimental point of view, time duration is 
not a problem given that, for example, a measurement 
that lasts 5 seconds can be easily performed. However, 
in the numerical framework the CFL stability condition 
pose severe restrictions on the time step At to be used, so 
that for intricate and large computational domains, a 5 
seconds event may involve several hours of computational 
time. 

Consider for example, the radiation impedance com- 
putation of the circular /a/ in two cases: i) a lossless 
impedance duct with /i = and ii) a lossy impedance 
duct with ji = 0.01. Let us first focus on the acoustic 
pressure collected at the first virtual microphone #1 (see 
Fig. 3) and plot its time evolution for the lossless and 
lossy cases in Fig. 4. For the former, the 30 ms duration 
of the simulation has not sufficed to attenuate the signal 
acoustic pressure inside the impedance duct, whilst it has 
decayed in about 15 ms for the lossy case. It should be 
remarked that given that the acoustic pressure has to be 
Fourier transformed to apply the TMTF method, it is 
necessary for it to vanish to zero during the simulation 
interval, to avoid spurious errors in this operation. Some 
more insight on what is going on inside the impedance 
duct can be obtained from the spectrograms (time vs 
frequency) of the acoustic pressure signals for the loss- 
less and lossy cases (see Fig. 5). As observed, the prob- 
lem arises from the difficulty to attenuate the first duct 
eigenmode. This is attributed to the fact that radiation 
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FIG. 4. Acoustic pressure evolution for circular /a/ at vir- 
tual microphone #1 for the lossless (jj>=0) and lossy (/x=0.01) 



is by far a more effective energy dissipating mechanism 
at high frequencies than at low frequencies. Therefore, 
the inclusion of wall damping clearly helps to overcome 
this problem and noticeably shortens the duration of the 
simulation. 



B. Analytical expressions for the complex axial wavenumber 

Introducing artificial damping at the impedance duct 
walls of the numerical model involves dealing with com- 
plex axial wavenumbers k z in expression (2) of the TMTF 
method. In experimental measurements, a calibration 
procedure is usually conducted for estimating the atten- 
uation factor (e.g. Boonen et ai, 2009). In the numer- 
ical case, analytical expressions can be deduced for k z . 
We will next show how to relate k z with the admit- 
tance boundary coefficient \i z in problem (4), for three- 
dimensional ducts of elliptical cross section. 

To proceed, it is convenient to express the wave 
equation (read also the Helmholtz equation) in ellip- 
tic cylindrical coordinates (£,77,2) (see e.g., Lowson and 
Baskaran, 1975). For each constant value of z, the coordi- 
nate lines correspond to confocal ellipses and hyperbolae. 
Curves of constant £ are ellipses whilst curves of constant 
r] are hyperbolae. Using this system of coordinates, the 
boundary condition (4d) at the duct wall (£ = £0) can be 
expressed as (see e.g., Oliveira and Gil, 2010) 



dp 



e 2 cos 2 r] p ■ 



at£ = £ , (6) 



where e = y/l — h\ja\ is the eccentricity of the el- 
lipse defining the duct boundary at £0, which has fo- 
cal distance f = a e e and major and minor semi-axes 
a e = /cosh£o and b e = /sinh£o. ft = k^a e is the re- 
duced (adimensional) wavenumber. To simplify notation, 
use will be also made of the parameter 



k±f 



k±a e e 



(7) 
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with k± = y/k^ — kl standing for the transverse 
wavenumber. 

Separation of variables for the Helmholtz equation in 
elliptic cylindrical coordinates results in the so-called 
Mathieu radial and angular equations. The solutions 
need to be 2tt periodic in 77, which plays the role some- 
how analogous to for the circular case (developed e.g. 
in Munjal, 1987). The periodic solutions are given by 
products of cosine elliptic functions ce m (77,(7), which are 
even, with the radial Mathieu functions Je m (£,q) re- 
lated to them, and by sine elliptic functions se m (ry, q), 
which are odd, and their corresponding radial Mathieu 
functions Jo m (£,q) (see Gutierrez- Vega, 2000). How- 
ever, no linear combination of ce m (r],q)Je m (£,q) and 
se m (r),q)Jo m (£,q) is allowed since the sets of character- 
istic values for ce m (77,(7) and se m (77,(7) are different (the 
situation is somehow similar to what occurs for rectan- 
gular ducts, see e.g., Munjal, 1987). Given that we are 
interested in the case of plane wave propagation inside 
the duct, we will be interested in the lowest modal in- 
dices and have to deal with the even case. 

The even radial Mathieu functions Je m (£,q) admit a 
series factorization in terms of Bessel functions. The fac- 
torization depends on m being even or odd so we can 
get even-even Je2&(£, q) and even-odd Je2/c+i(£, q) radial 
Mathieu function developments. As we want to consider 
the case m = we have to deal with the former. It turns 
that (Gutierrez- Vega, 2000) 



ce 2 k(0,q) 

An 



J2A 2j J2 j (2 y /qsmhO (8) 



3=0 



with J m corresponding to the Bessel function of order m. 
Taking k = j = we get 

p(€,rj,z) ~ ce (r],q)ceo(0,q)Jo (2y/qswh£) , (9) 

with derivative 

dp(£,ri,z) 



-ce (f],q)ce (0,q) 
x 2 v /</cosh£Ji (2y/qswh€) . (10) 



Substituting (9) and (10) into the boundary condition 
(6) and evaluating at £ = £0 we get 

2^cosh^ ( 2 ^ si nhCo) ^— 
ft J (2^ sum £ J 

In order to eliminate the 77 dependence in the r.h.s of (11) 
we can integrate at both sides from to tt/2. Changing 
variables to v = cos 77 shows that the integral in the r.h.s 
is a complete elliptic integral of the second kind, whose 
solution can be expressed as a series in terms of even 
powers of the eccentricity (see Abramowitz and Stegun, 
1970) 



r/ 2 r 

/ v 1 — e 2 cos 2 rjdr] = / 
Jo Jo 



1 VT 



00 
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-dv 



/(e) (12) 



Using this result, approximating the Bessel functions to 
first order, Jo(x) ~ 1, J\(x) ~ x/2, and considering the 
notation introduced above, it follows from Eq. (11) that 



jfjL z Kl(e) 



.\± z knl{e)a\e 2 



7r cosh £0 sinh £0 
From (7) we then get 

AVzkoI(e) 



irb e 



7rb e 



(13) 



(14) 



and the axial wavenumber that we were looking for be- 
comes 



,W(e) 



(15) 



In the case of zero eccentricity we get 1(0) = 7r/2, a e = 
b e = a, and recover the axial wavenumber for the circular 

case 



(16) 
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FIG. 6. Radiation resistance R r and reactance X r for vowel /i/ with circular (top) and elliptical (bottom) mouth aperture 
for different pairs of virtual microphones. (a)(c) Real wavenumber computations. (b)(d) Complex wavenumber computations. 
Solid lines correspond to FEM simulations and square symbols to the theoretical model, only available for a circular piston. 



which can be derived from first order approximations in 
the characteristic equation for circular cylindrical ducts 
in Munjal (1987). 

On the other hand, note that to compute the 
impedance of vocal tracts, impedance ducts of differ- 
ent cross sections and shapes will be necessary. This 
implies that the total area of the impedance duct walls 
will be different in each case, so that if we use the same 
boundary admittance coefficient \i z we will get different 
internal dissipation. To guarantee the same amount of 
absorption, we must ensure that the axial wavenumber 
remains the same for all cases. Suppose that we want to 
get the same dissipation in an elliptical impedance duct i 
than in another elliptical impedance duct j, with respec- 
tive eccentricities ej, and semi-minor axes 6 e ^, b e j. 
Equating their axial wavenumber s (15) we can compute 
fi Zi i from /jL z j as 

As mentioned in section III.A, in this work a boundary 
admittance coefficient of \i z = 0.01 has been chosen for 
the circular vowel /a/. Therefore, to introduce the same 
amount of dissipation in all simulations, we have made 
use of (17), with j =/a/ and i =/i/,/u/. 



C. Accuracy of complex axial wavenumbers 

Several assumptions have been made in order to derive 
the axial complex wavenumbers in (15) and (16). Let us 
next check if these simple expressions are precise enough 
for vocal tract impedance computations and what would 
happen if real wavenumbers were considered instead. To 
do so, we have computed the radiation impedance of 
the circular and elliptical /i/ using complex and real 
wavenumbers, and considered different pairs of virtual 
microphones, located at different distances from the duct 
exit. The computed radiation impedances, splitted in 
terms of resistance R r and reactance X r , have been plot- 
ted in Fig. 6 and compared to the theoretical model of the 
baffle set in a spherical surface (Morse and Ingard, 1968), 
for the circular case. The single point arrow on the R r 
curves indicates their tendency when the selected pairs of 
microphones used for the computations are moved away 
from the duct exit (see duct scheme in the upper left side 
of figures). The double point arrow indicates that the X r 
curves do not follow any clear tendency when changing 
the pairs of virtual microphones. 

If we have a look at Eq. (2), we observe that the 
term e^ 2kzXl is used to translate the reflection coef- 
ficient R\ to the reference plane (x\ is the distance 
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from the first virtual microphone to the reference plane). 
Taking into account the first order approximation (1 + 
x) 1 / 2 ~ 1 + x/2 in the expression for the complex 
wavenumber in (15) (which is valid for hard wall be- 
havior \Z W \ ^> poco), it follows that fc^oo = k z ~ 
ko — j2ji z I{e) / (jibe). Substituting into the propagator 

e j2k zXl results in e j2x 1 ko e j2x 1 ^ z I(e)/( 7 rb e )^ SQ ^ the fac _ 

tor e J 2x iVzi(e)/(7rb e ) ]j e m i S sing if a real wavenumber 
is used. Note that the higher the value of x\ the larger 
will be the error. This can be appreciated in Fig. 6a 
for the circular case, where R r clearly departs from the 
theoretical curve as x\ increases. As opposite, the error 
is inexistent when the appropriate complex wavenumber 
of equation (16) is used in (2), and the computed curve 
perfectly matches the theoretical one (see Fig. 6b). The 
reactance X r seems not to be as much affected as the 
resistance R r by the error. 

A very similar behavior can be observed for the ellip- 
tical case as shown in Fig. 6c and Fig. 6d. Besides, note 
that both for the circular and elliptical cases, some dis- 
crepancies are found for the reactance values at the high 
frequency range, even when using the correct complex 
wavenumber s. This is attributed to numerical errors in 
the simulations to be discussed in section III.E. 



D. Plane wave propagation restriction 

One of the main frequency limitations of the TMTF 
method is that of plane wave propagation assumption. 
Next we will examine how this condition can be relaxed 
for vocal tract impedance ducts with circular and ellipti- 
cal sections. 

For a circular impedance duct of radius a, the first 
non-planar eigenmode is the (1,0) mode, with a cut- 
off frequency /(i,o) = 1.84co/(27ra) (Fletcher and Ross- 
ing, 1988). The next one, is the (2,0) mode, with 
/(2,o) = 3.05co/(27ra), followed by the (0,1) mode, with 
/(o,i) = 3.80co/(27ra). In the experimental framework, 
the (1,0) mode limits the working frequency range of the 
TMTF method. However, in the numerical framework, 
this limitation can be overcome by a proper location of 
the virtual microphones. If we have a look at the pressure 
distribution of the first three eigenmodes (see Fig. 7a), 
we can observe nodal planes (corresponding to lines in 
the duct sections of the figure) that cross the centerline 
of the duct for the (1,0) and (2,0) modes. Consequently, 
if we precisely locate the virtual microphones at the cen- 
terline, the pressure there will not be affected by these 
modes, but will be only due to plane wave propagation. 
The TMTF method will be then still applicable at these 
frequencies and limited by the (0,1) mode, which is the 
first eigenmode that does not have a nodal line at the 
center of the duct (see Fig. 7a). 

With regards to elliptical ducts, the frequency range 
can also be extended. In Fig. 7b we present the pressure 
patterns of the first four eigenmodes for the elliptical /a/ 
impedance duct (these can be computed e.g., from the 
formulas in Oliveira and Gil, 2010). The first non-planar 
mode is the even mode (1,1). This mode has a nodal 
plane in the center, so that locating again the virtual 




(a) (b) 



FIG. 7. (Color online) Plots of the lower standing modes in 
an impedance duct with (a) circular and (b) elliptical section. 
Discontinuous lines represent nodal planes. 



microphones in the centerline of the impedance duct will 
allow to extend the analysis up to the frequency of the 
even mode (2,1), which is the limiting one in this case. 

TABLE I. Frequency values in KHz for the first eigenmodes of 
the circular and elliptical impedance ducts used to compute 
the radiation impedance Z r and the input impedance Zi n for 
vowels /a/, /i/ and /u/. Stars denote maximum frequency 
of analysis when the virtual microphones are located at the 
centerline of the impedance duct. 



Circular 




Elliptical 




(0,1) (1,0) (1,1) 


E(l,l) 


E(2,l) 


E(3,l) 


0(1,1) 


/a/ 8.21 13.62 17* 


4.68 


8.6* 


12.47 


13.79 


Z r /i/ 18.8 31.13 38.8* 


8.27 


15.27* 


22.18 


39.3 


/u/ 44.77 74.2 92.46* 


22.23 


40.98* 


59.49 


84.2 


Z in /a/ 23.72 39.32 48.99* 


13.57 


24.97* 


36.2 


40.03 



In Table I we have computed the first modes of the 
impedance ducts used in this work, with stars denot- 
ing in each case the first mode without a nodal plane 
at its center (i.e., the limiting one). As expected, the 
impedance duct with a strongest restriction is that of 
vowel /a/ (largest mouth aperture), whereas vowel /u/ 
(smallest mouth aperture) presents the less stringent con- 
dition. It can also be observed that working with ellipti- 
cal mouth apertures results in more restrictive frequency 
ranges. Note however that by locating the virtual mi- 
crophones at the centerline of the impedance duct, the 
upper frequency limit of the experimental TMTF method 
has been almost increased by a factor ~ 2. Except for 
the elliptical /a/, the desired frequency range of analysis 
that goes up to fmax — 10 KHz can thus be attained. 



E. The singularity of the TMTF method 

It is well-known that a singularity occurs in the experi- 
mental TMTF method when half the wavelength of the 
acoustic pressure equals the microphone spacing s, or 
one of its multiples. The minimum frequency that satis- 
fies this condition is termed the critical frequency, f cr = 
Co/ (2s), and imposes a high frequency limit f u < f cr to 
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Frequency [KHz] 

FIG. 8. Resistance and reactance of the elliptical /i/ for dif- 
ferent microphone spacings s. The dashed line stands for 
s = 0.1 cm (equal to the mesh size h), the dotted line for 
s = 1.5 cm (close to the singularity), and solid lines corre- 
spond to intermediate values of s. 



the method. On the other hand, a minimum low fre- 
quency limit fd also exists, because the pressure differ- 
ences measured by the two microphones will be negli- 
gible if their distance apart is very small compared to 
the measured wavelength. Working close to both lim- 
its is not recommended and a suitable option is to take 
f d = 0.2f cr < f < O.Sf cr = f u (as proposed in Jang 
and Ih, 1998). 

The above situation contrasts with that of numerical 
simulations. Suppose that we are interested in making 
a simulation up to a given frequency fmax with corre- 
sponding wavelength \ m in- If we had unlimited com- 
putational resources, so that we could work with a fine 
enough mesh to meet the standard criterion of 10 nodes 
per wavelength for Xmin (see e.g., Ihlenburg, 1998), nu- 
merical errors would be negligible. In such situation the 
only frequency upper limitation would come from the 
plane wave restriction discussed in section III.D, and we 
could take e.g., fd = 0<f< fmax = for = fu as the 
operational working frequency range. This expression to- 
gether with f cr = Co/ (2s) provides the following possible 
values for the virtual microphone spacing s 

h<s< 0.5A min . (18) 

Note from (18), that obviously the microphone spacing s 
cannot be smaller than the mesh size h. 

However, working with very fine meshes can be unfeasi- 
ble in many simulations and the 10 nodes per wavelength 
criterion is often sacrificed to lower the computational 
cost. This results in the appearance of some numerical 
errors at high frequencies as happens with the simulations 
performed throughout this work, where 34.5, 6.9 and 4.6 
nodes per X m in have been respectively taken within the 
impedance duct, the outer volume and the PML region 
(see section II. D for mesh details). In Fig. 8 we have plot- 
ted the radiation resistance and reactance for the ellip- 
tical jij computed using different microphone spacings, 
s = 0.1, 0.25, 0.5, 0.75, 1, 1.25, 1.5 cm, which respectively 
yield s/Xmin ~ 0.03, 0.07, 0.15, 0.22, 0.29, 0.36, 0.43. 



0.1 0.2 0.3 0.4 0.5 

s/X 

FIG. 9. Singularity factor SF of the TMTF method for a 
lossless impedance duct. The circle and the square denote 
the low and high limits for the microphone spacing s obtained 
when SF = 1.7. Vertical lines correspond to the microphone 
spacing configurations used in the example of Fig. 8. 



The first spacing (s = 0.1 cm) has been chosen equal to 
the mesh size h inside the impedance duct (dashed curve 
in Fig. 8), whilst the last one (s = 1.5 cm) has a value 
close to the singularity value 0.5A m i n (dotted curve in 
Fig. 8). No differences can be appreciated in Fig. 8 be- 
tween the various R r and X r curves for the low-mid fre- 
quency range as there are no significant numerical errors 
at these frequencies for the present simulations. Slight 
differences only become apparent for frequencies higher 
than 8 KHz (see zoom in Fig. 8), the largest ones pre- 
cisely corresponding to the reactances computed with the 
limiting values of s = 0.1 cm (mesh size) and s = 1.5 cm 
(singularity value). 

The above example shows that the condition in 
Eq. (18) for the microphone spacing is too loose and that 
stronger requirements are needed in practice for the spac- 
ing limiting values. In order to define them we will re- 
sort to the so called singularity factor (SF) of the TMTF 
method, introduced by Jang and Ih (1998). The SF in- 
dicates the sensitivity of the TMTF method to errors in 
the input pressures P\(f) and ^(Z); the higher the SF 
value the stronger the influence of the error sources in 
the computed impedances. When computing the SF in 
the experimental framework, it is assumed that all errors 
are of the white type, uncorrelated and with constant 
variance. This will not be the case of a single numeri- 
cal simulation, the error being totally deterministic i.e., 
always the same if we repeat the computation. How- 
ever, if we consider the simulation of a given vocal tract 
impedance being representative of a certain ensemble av- 
erage of vocal tracts having e.g., slight different geometry 
details and material characteristics, we can hypothesize 
that errors arising from these simulations would satisfy 
the error requirements for the SF computation. It would 
then be logical to demand that the chosen spacing for the 
virtual microphones results in a small SF value. 
In Fig. 9, we have plotted the SF curve (thick line) for 
the standard TMTF method in a lossless impedance duct 
(implementation of the SF for a lossy impedance duct is 
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out of the scope of this work), according to the proce- 
dure described in Jang and Ih (1998). However, instead 
of representing the SF dependence with frequency, we 
have plotted it against the microphone spacing s (nor- 
malized by X m i n ). Moreover, we have also plotted the 
microphone spacings corresponding to the resistance and 
reactance curves in Fig. 8 as vertical lines. As observed 
in the figure, the SF values for the extreme spacings cor- 
responding to the mesh size and singularity values are 
much higher than the threshold value of SF < 1.7, rec- 
ommended in the experimental framework (see Jang and 
Ih, 1998). According to this criterion we can get a more 
restrictive range for the virtual microphone spacing than 
that provided by Eq. (18), namely 

ti < 0.1A min < s < 0AX min . (19) 

The optimum, and thus recommended, spacing s will be 
that minimizing SF which is close to s ^ 0.25A m i n in 
Fig. 9. In our computations, we have used s = 1 cm (see 
section II. D) that yields s ~ 0.29A min . 

IV. VOCAL TRACT IMPEDANCES 

In this section, radiation and input impedances for vowel 
vocal tracts using the described FEM-TMTF approach 
will be presented. As detailed in section II, we have 
computed the radiation impedance for vowels /a/, /i/ 
and /u/, considering circular and elliptical mouth aper- 
tures, and the input impedance of vowel /a/ for circular 
and elliptical vocal tracts. 

In Fig. 10 we present the results for the radiation 
impedance, split in terms of the radiation resistance 
(Fig. 10a) and reactance (Fig. 10b). All results are pro- 
vided up to 10 KHz except for the elliptical /a/, the 
analysis only being valid in this case up to ~ 8 KHz 
because of the limiting even mode (2,1) (/ = 8.6 KHz, 
see Table I). Note however, that we can reach this value 
thanks to proper location of the virtual microphones (see 
section III.D) and that the experimental TMTF would 
have only let us to measure impedance values up to 
/ = 4.2 KHz, corresponding to the even mode (0,1) (see 



Table I). Similarly, the plane wave frequency restrictions 
for the circular /a/ (/ = 8.21 KHz) and for the ellipti- 
cal /i/ (f = 8.27 KHz) in the experimental TMTF, can 
be easily overcome thanks to centerline microphone po- 
sitioning. Thus, for all analyzed cases, except for the 
elliptical /a/, the radiation impedance can be computed 
for the whole frequency range of interest (0 10] KHz, 
without problems (see Table I). 

The differences between the resistance and reactance 
curves for the circular and elliptical mouth apertures in 
Fig. 10 can be justified as follows (remember that for a 
given vowel only the shape of the mouth, elliptical or 
circular, changes but not the total mouth surface). In 
the plane wave propagation regime, the curves should be 
almost identical. However, as explained in section II. D, 
the reference surface for elliptical mouth apertures is de- 
fined at the intersection between the major semi- axis of 
the impedance duct with the sphere. This implies that 
there is small portion of the duct beyond the reference 
surface, which reaches the intersection between the minor 
semi- axis and the sphere, and plays somehow the role of 
some kind of lips. Its effects on the radiation impedance 
of the elliptical ducts can be understood as those of a 
modified end correction. To first order approximation, 
when R% + X% « 1, the inertial end correction only 
affects the reactance. This is what can be observed for 
vowel /u/, which shows the same resistance values for the 
circular and elliptical cases (see Fig. 10a), but different 
ones for the reactance (see Fig. 10b). For vowels /a/ and 
the condition B% + X% << 1 is no longer satisfied, 
say for frequencies bigger than 4 KHz, and differences in 
resistance values can become more clearly appreciated. 
It is to be noted that above 4 KHz, differences between 
the elliptical and circular cases are not only due to the 
end correction effects but also to the presence of higher 
order modes in the impedance ducts. 

Finally, we have computed the input impedance of the 
circular and elliptical /a/ and plot their moduli in Fig. 11. 
In contrast to the radiation impedance computation, the 
input impedance for elliptical /a/ can be computed with 
no problem up to 10 KHz. This is so because the glot- 
tal section is much smaller than the mouth aperture; a 
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FIG. 11. Input impedance for the circular and elliptical /a/. 



narrower impedance duct is then required, which will 
have the first non-planar mode beyond 10 KHz. Look- 
ing at Fig. 11, we observe that below 5 KHz the input 
impedances moduli of the circular and elliptical /a/ are 
very similar indicating that we are in the plane wave 
propagation range. Note however, that there is a cer- 
tain formant shift to lower frequencies for the elliptical 
/a/, because its radiation reactance is higher than that 
of the circular /a/ (see Fig. 10b). Above 5 KHz, dif- 
ferences become notorious due to non-planar high order 
mode effects. 



V. CONCLUSIONS 

In this paper, we have shown how the experimental two- 
microphone transfer function method (TMTF) can be 
adapted to compute impedances of elliptical vocal tracts, 
from time domain simulations of the irreducible wave 
equation. This avoids having to compute the acoustic 
velocity field to do so. 

However, using the TMTF in the numerical context 
it is not straightforward and several considerations have 
been analyzed. These have resulted in the following ob- 
servations. First, it is important to impose losses at the 
impedance duct walls to attenuate the first duct eigen- 
mode and achieve reasonable computational times. This 
implies using complex wavenumbers in the TMTF ex- 
pressions, which have been derived for three-dimensional 
cylindrical impedance ducts with elliptical cross section. 
Second, we have seen that the frequency range of va- 
lidity of the experimental TMTF method can be al- 
most doubled by locating the virtual microphones at the 
impedance duct center line. This allows overcoming to 
some extent the plane wave limitation of the experimen- 
tal TMTF. Third, a range of possible values for the vir- 
tual microphone spacing has been proposed making use 
of the so called singularity factor. An optimum value of 
a quarter wavelength of the maximum frequency to be 
solved is recommended. 

The radiation impedance of vowels /a/, /i/ and /u/ 
with circular and elliptical mouth apertures, and the in- 
put impedance of vowel /a/, again for the elliptical and 



circular cases have been computed to test the proposed 
approach. Simple geometries of vocal tracts available in 
the literature have been used though any realistic one 
(e.g., generated from MRI) could have been used instead. 
All simulations have been carried out using a finite ele- 
ment (FEM) approach with a perfectly matched layer 
(PML) to allow outward waves from the mouth propa- 
gate to infinity. 
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APPENDIX A: TIME DOMAIN FEM APPROACH 

In this appendix we include the details of the time domain 
FEM approach used to solve problem (4) in section I. 
This consists in extending the finite difference PML for- 
mulation for the irreducible wave equation in Grote and 
Sim (2010), to the FEM framework. 

If we replace the Sommerfeld radiation condition with 
the PML, the original acoustic wave equation (4) be- 
comes modified to 



d t <t>i 
d t tp - 



- V • <fi — adtp — Pp — 7^ 
claidip + clbidiip, Vi 



(Ala) 



1,2,3 (Alb) 
(Ale) 



in ft, t > 0, with boundary and initial conditions 

Vp • n = -po/Sd t Q = g on T G , t > 0, (Aid) 

Vp • n = —fi/codtp on Tw, t > 0, (Ale) 

Vp n = on T H , t > 0, (Alf) 

Vp • n = on r^, t > 0, (Alg) 

p = 0, d t p = in ft, t = 0, (Alh) 

(j>i = 0, d t <t>i = 0, Vi = 1, 2,3 in ft, t = 0, (Ali) 

^ = 0, d t ^ = inft,£ = 0. (Alj) 

Note that (Ala) modifies the right hand side (r.h.s.) of 
(4a) with the inclusion of some extra terms involving the 
auxiliary functions ip and = (0i, 02,^3). Moreover, 
four additional scalar equations (Alb)-(Alc) are needed 
for these additional functions. In what concerns the co- 
efficients a, /?, 7, a^, bi in (Al), they depend on the 
damping profiles & and are given by a = £i + £2 + £3, 

P £16 + 66 + 66, 7 = 666, 0*1 =6 + 6-6, 
q>2 6+6-6, «3 6+6-6, h = 66, h 66 and 

bs = 66- The damping profiles 6 are used to control 
the amount of absorption in the PML and many options 
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exist for them. Following Grote and Sim (2010), we have 
used 



■ { 2ir\xi-li\\ 



2tt 



(A2) 



for li < \xi\ < h+ Li. & is a constant accounting for the 
damping effect in the i-th direction, k is the i-th coordi- 
nate of the PML layer and Li the thickness of the PML 
region in the i-th direction. The constant ^ depends on 
the discretization and thickness of the layer and can be 
computed as 



log 



(A3) 



with standing for the relative reflection at the bound- 
ary of the PML. The PML boundary can be truncated 
using either a Dirichlet or a Neumann homogeneous con- 
dition (the latter has been our choice, see (Alg)). On 
the other hand, notice that for li < \xi\, = 0, i.e., 

outside the PML, the damping profiles & vanish. In this 
case, the modified PML wave equation (Ala) reduces to 
the standard wave equation in (4a). 

The set of partial differential equations (Ala)- (Ale) 
supplemented with boundary and initial conditions 
(Ald)-(Alj) constitute the problem we aim at solving. 
If a FEM approach is to be used to find a numerical 
solution for (Al), we first have to set this problem in 
its weak or variational form. As usual, we first multi- 
ply (Al) by test functions q,Vi,w (q for the pressure, Vi 
for the first auxiliary functions 0^, and w for the second 
auxiliary function ip) and we integrate over the computa- 
tional domain Q. Applying the divergence theorem and 
making use of boundary conditions, we get the weak form 
of the problem, which consists in finding p, fa and ip such 
that 

(<3S dip) + c (q, ^d t p) Fw + c\ (V<?, Vp) = c\ (q, g) Fc 

3 

i=\ 

(A4a) 



- (q,Pp) - (<7,7^), 

(vi, d t (j>i) - (vi,£i<t>i) + Cq (vi, didip) 
+ c 2 (v i: bidi^), Vi = l,2,3, 
(w,dti/>) = (w,p) , 

in 17, t > 0, with initial conditions 



(A4b) 
(A4c) 



(q,p) = 0, (q,d t p) = 0, (A4d) 
(v^ fa) = 0, (v u d t fa) = 0, Vi = 1, 2, 3, (A4e) 
K^)=0, (w,dt*l>) = 0, (A4f) 

in t = 0, for all v^, w. To shorten notation, in (A4) 
the integral of the product of any two functions f,g'mQ 
has been written as 



(/,<?) 



fgdfl, 



(A5) 



whilst integrals over boundaries have been explicitly in- 
dicated e.g., (f,g) FG . 

To find a numerical solution to problem (A4), we have 
to discretize it both in space and time. Let us first 
discretize it in space using a FEM formulation. Given 
a finite element partition of £1 with n e \ elements and 
rip nodes and discretizing problem (A4) following the 
Galerkin method, we have expanded the discrete ver- 
sions of the unknowns p, fa, ip and the discrete ver- 
sions of the test functions q, vi and w in piecewise lin- 
ear terms of piecewise linear shape functions N(x) (e.g., 
Ph = Yl^Li ^ h P h ) ^ so that we get the following time 
evolving algebraic matrix system 

MP + c BP + clKP = c 2 L 

3 

+ ~ ~ M P P ~ M 7 *, (A6a) 

i=i 

M<$>i = -Mfr&i + ^Bi^P 
+ c 2 B iM V, Vi = 1,2,3, 

4/ = p, 

in t > 0, with initial conditions 

P = 0, P = 0, 
$i =0, *< = 0, Mi = 
^ = 0, ^ = 0, 



1,2,3, 



(A6b) 
(A6c) 

(A6d) 
(A6e) 
(A6f) 



in Q,t = 0. In (A6), P, &i and \l/ stand for the vectors 
of nodal values that respectively correspond to the pres- 
sure and auxiliary functions (e.g., P = (P 1 • • • P Up ) T ), 
whereas the remaining matrix and vector entries are 



iven by 










M ab = 


(N a ,N b ) , 


M ab = 


(N a ,a h N b ), 


(A7a) 


Mf = 


(N a ,/3 h N b ), 


M° b = 


(N a , lh N b ), 


(A7b) 


Mf = 

si 


(N a ,UN b ), 


B ab = 


(N*,»N b ) rw , 


(A7c) 


nab 

B i,ai - 


{N a ,a th d t N b ), 


Bf = 


(N a ,m b ), 


(A7d) 


Tjab 

U iM = 


{N a ,b ih diN b ), 


R ab = 


(VN a ,VN b ) , 


(A7e) 


L a = 


(N a ,g) rc . 






(A7f) 



As usual, the domain integrals in (A7) are to be under- 
stood as the summation of integrals over elements Q e , 
i.e., = Ee=i( , ' , ) fi e- The subscripts h denote the 

discrete versions of the corresponding continuous vari- 
ables. 

Let us next proceed to the time discretization of (A6). 
A finite difference approach has been used to do so. We 
consider a constant time step (At = t n+1 — t n ) discretiza- 
tion of the time interval (0, T) into < t 1 < • • • < t 



.71—1 



< 



f 1 < t 



n+1 



< 



< t 



N 



T. A second order finite dif- 



ference central scheme has then been implemented for 
the pressure time derivatives, whilst a first order central 
scheme has been used for the time derivatives of the aux- 
iliary variables in the PML region. This is so because first 
order schemes are known to introduce strong numerical 
dissipation, which in this case is advantageous to help 
absorbing the waves crossing the PML. 
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Inserting these time derivative approximations in (A6) 
yields 

pn+l 2P n + P n ~^ pn+1 pn—1 

M — + cqB- 



At 2 



2At 



c\KP n = c 2 L n + Bi$? - M a 



pn-\-l pn—1 



2At 



- M R P n - M, 



(A8a) 



M- 



At 



= -M P 



n+l 



<I> r 



clB. 



>n+l 



pn 



•0- £ -%ai 



n+1/2 



^,n+l/2 _ q>n-l/2 



At 



= P™, 



(A8b) 



(A8c) 



with initial conditions 



p° 


= 0, P 1 


= 0, 


(A8d) 


I 


= o, *! 


= 0, Mi = 1,2,3, 


(A8e) 


^0 


= 0, tf 1 


= 0. 


(A8f) 



Note that (A8) corresponds to a purely explicit scheme, 
where all the unknowns can be calculated at time step 
n + l from values already known from previous steps. 
Solving (A8) at time step t = n + 1 involves 

1. Use (A8c) to compute ^ n+1 /2 

2. Insert ^ n+1 /2 mto (A8a) and compute P n+1 

3. Update using (A8b), ^+V2 and pn+i 

To fasten all computations, matrix inversion is avoided as 
usual by means of a lumped approximation for all mass 
matrices (A7a)-(A7c) (Hughes, 2000). 

The system of equations (A8) constitutes the final nu- 
merical scheme that has been used for all the computa- 
tional simulations in this work. 
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